R Notebook

Load packages in

pacman::p_load(sf, readr, tidyverse, tmap)

Load Data in

Geospatial Data

seoul <- st_read(dsn = "data_final/seoul_adm2",
                 layer= "seoul_adm2") %>%
  st_transform(crs = 5179)

Reading layer `seoul_adm2' from data source 
  `/Users/junhaoteo/Documents/junhao2309/Geospatial-Analytics-Project-KR/data_final/seoul_adm2' using driver `ESRI Shapefile'
Simple feature collection with 748 features and 6 fields
Geometry type: MULTIPOLYGON
Dimension:     XY
Bounding box:  xmin: 126.7644 ymin: 37.42854 xmax: 127.1838 ymax: 37.70139
Geodetic CRS:  WGS 84

clinics <- st_read(dsn = "data_final/seoul_clinics",
                   layer = "seoul_clinics") %>%
  st_transform(crs = 5179)

Reading layer `seoul_clinics' from data source 
  `/Users/junhaoteo/Documents/junhao2309/Geospatial-Analytics-Project-KR/data_final/seoul_clinics' using driver `ESRI Shapefile'
Simple feature collection with 40 features and 4 fields
Geometry type: MULTIPOINT
Dimension:     XY
Bounding box:  xmin: 126.8894 ymin: 37.48397 xmax: 127.0858 ymax: 37.63072
Geodetic CRS:  WGS 84

hospitals <- st_read(dsn = "data_final/seoul_hospitals",
                     layer = "seoul_hospitals") %>%
  st_transform(crs = 5179)

Reading layer `seoul_hospitals' from data source 
  `/Users/junhaoteo/Documents/junhao2309/Geospatial-Analytics-Project-KR/data_final/seoul_hospitals' using driver `ESRI Shapefile'
Simple feature collection with 88 features and 4 fields
Geometry type: MULTIPOLYGON
Dimension:     XY
Bounding box:  xmin: 126.8351 ymin: 37.45503 xmax: 127.1581 ymax: 37.64912
Geodetic CRS:  WGS 84

residence <- st_read(dsn = "data_final/seoul_residential",
                     layer = "seoul_residential") %>%
  st_transform(crs = 5179)

Reading layer `seoul_residential' from data source 
  `/Users/junhaoteo/Documents/junhao2309/Geospatial-Analytics-Project-KR/data_final/seoul_residential' using driver `ESRI Shapefile'
Simple feature collection with 9640 features and 5 fields
Geometry type: MULTIPOLYGON
Dimension:     XY
Bounding box:  xmin: 126.8081 ymin: 37.43452 xmax: 127.1722 ymax: 37.67945
Geodetic CRS:  WGS 84

road <- st_read(dsn="data_final/seoul_roads",
                layer = "seoul_roads_v2")%>%
  st_transform(crs = 5179)

Reading layer `seoul_roads_v2' from data source 
  `/Users/junhaoteo/Documents/junhao2309/Geospatial-Analytics-Project-KR/data_final/seoul_roads' using driver `ESRI Shapefile'
Simple feature collection with 100303 features and 10 fields
Geometry type: MULTILINESTRING
Dimension:     XY
Bounding box:  xmin: 126.7662 ymin: 37.42867 xmax: 127.1831 ymax: 37.70137
Geodetic CRS:  WGS 84

Select wanted datasets

clinic <- clinics %>%
  select(1,3)
hospitals <- hospitals %>%
  select(1,3)
residence <- residence %>%
  select(1,3)
road <- road %>%
  select(1,3)

Aspatial Data

population <- read_csv("data_final/population_seoul.csv")

Rows: 25 Columns: 21── Column specification ──────────────────────────────────────────────────────────────────────────────────────────────────────────
Delimiter: ","
chr  (1): By administrative divisions(eup, myeon, dong)
dbl (20): Total population (Person), Male (Person), Female (Person), Korean - total (Person), Korean - male (Person), Korean -...
ℹ Use `spec()` to retrieve the full column specification for this data.
ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.

Data Wrangling

Geospatial

pop <- population %>%
  select(1,2) %>%
  mutate(Elderly_pop = round(0.1758*`Total population (Person)`,0)) %>%
  rename(name_en =`By administrative divisions(eup, myeon, dong)`)

Gu <- pop$name_en

seoul <- seoul %>%
  filter(name_en %in% Gu) %>%
  select(name_en)

seoul_pop <- seoul %>%
  left_join(pop,
             by = "name_en") 

st_write(seoul_pop,
         dsn = "data_final",
         layer = "seoul_pop",
         driver = "ESRI Shapefile")

seoul_pop <- st_read(dsn = "data_final/",
                     layer = "seoul_pop")

tm_shape(seoul_pop)+
  tm_fill("Elderly_pop", 
          style = "quantile", 
          palette = "Greens",
          title = "Elderly Population") +
  tm_layout(main.title = "Elderly Population",
            main.title.position = "center",
            main.title.size = 1.2,
            legend.height = 0.45, 
            legend.width = 0.35,
            frame = TRUE) +
  tm_borders(alpha = 0.5) +
  tm_compass(type="8star", size = 2) +
  tm_scale_bar() +
  tm_grid(alpha =0.2) +
  tm_credits("Source: Data.korea", 
             position = c("left", "bottom")) +
  tmap_style("classic")

tmap style set to "classic"
other available styles are: "white", "gray", "natural", "cobalt", "col_blind", "albatross", "beaver", "bw", "watercolor" 
Credits not supported in view mode.
Compass not supported in view mode.

Start of Network Constrained Spatial Point

Packages needed

pacman::p_load(sf, spNetwork)

tmap_mode("view")

tmap mode set to interactive viewing

tm_shape(roads_by_gwanak) +
  tm_lines()

roads_by_gwanak <- roads_by_gu_new %>%
  filter(name_en == "Gwanak-gu")
hospital_by_gwanak <- hospital_by_gu %>%
  filter(name_en == "Gwanak-gu")
residence_by_gwanak <- residence_by_gu %>%
  filter(name_en == "Gwanak-gu")

Ensure datasets are in the right format for NKDE

hospital_by_gu <- st_centroid(hospital_by_gu)

Warning: st_centroid assumes attributes are constant over geometries of x

residence_by_gu <- st_centroid(residence_by_gu)

Warning: st_centroid assumes attributes are constant over geometries of x

if ("MULTILINESTRING" %in% st_geometry_type(roads_by_gu)){
  converted <- st_cast(roads_by_gu[which(st_geometry_type(roads_by_gu) == "MULTILINESTRING"),], "LINESTRING")
  linestring <- roads_by_gu[which(st_geometry_type(roads_by_gu) == "LINESTRING"),]
  roads_by_gu_new<- rbind(linestring, converted)
}

Warning: repeating attributes for all sub-geometries for which they may not be constant

roads_by_gwanak <- roads_by_gu_new %>%
  filter(name_en == "Gwanak-gu")
hospital_by_gwanak <- hospital_by_gu %>%
  filter(name_en == "Gwanak-gu")
residence_by_gwanak <- residence_by_gu %>%
  filter(name_en == "Gwanak-gu")

tmap_mode("view")
tm_shape(roads_by_gwanak) +
  tm_lines() +
  tm_shape(hospital_by_gwanak) +
  tm_dots(col = "green") +
  tm_shape(residence_by_gwanak) +
  tm_dots(col = "red")

Preparing lixels object

# Filter out only LINESTRING geometries
lixels <- lixelize_lines(roads_by_gwanak, 700, mindist = 350)

Generating line centre points

samples <- lines_center(lixels)

Performing NetKDE

hospital_by_gwanak_sp <- as_Spatial(hospital_by_gwanak)

densities <- nkde(roads_by_gwanak, 
                  events = hospital_by_gwanak,
                  w = rep(1,nrow(hospital_by_gwanak)),
                  samples = samples,
                  kernel_name = "quartic",
                  bw = 300, 
                  div= "bw", 
                  adaptive = FALSE,
                  method = "continuous", 
                  digits = 1, 
                  tol = 1,
                  grid_shape = c(1,1), 
                  max_depth = 8,
                  agg = 20, #we aggregate events within a 5m radius (faster calculation)
                  sparse = TRUE,
                  verbose = FALSE)

Visualising NETKDE

samples$density <- densities*1000
lixels$density <- densities*1000

tmap_mode('view')

tmap mode set to interactive viewing

tm_shape(lixels)+
  tm_lines(col="density")

Network constrained G and K function

Ho: The observed spatial point events (i.e distribution of childcare centres) are uniformly distributed over a street network in Punggol Planning Area.

kfun_hospital <- kfunctions(roads_by_gwanak, 
                             hospital_by_gwanak,
                             start = 0, 
                             end = 30, 
                             step = 1, 
                             width = 50, 
                             nsim = 50, 
                             resolution = 50,
                             verbose = FALSE, 
                             conf_int = 0.05)

kfun_hospital$plotk

Geographical accessibility

Using Euclidean Distance

library(sf)
library(rgeos)

distance_matrix <- as.matrix(st_distance(residence_by_gwanak, hospital_by_gwanak))